Ion specificity and anomalous electrokinetic effects in hydrophobic nanochannels 
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We demonstrate with computer simulations that anomalous electrokinetic effects, such as ion 
specificity and non-zero zeta potentials for uncharged surfaces, are generic features of electro-osmotic 
flow in hydrophobic channels. This behavior is due to the stronger attraction of larger ions to the 
"vapour-liquid-like" interface induced by a hydrophobic surface. An analytical model involving 
a modified Poisson-Boltzmann description for the ion density distributions is proposed, which al- 
lows the anomalous flow profiles to be predicted quantitatively. This description incorporates as 
a crucial component an ion-size-dependent hydrophobic solvation energy. These results provide an 
effective framework for predicting specific ion effects, with important implications for the modeling 
of biological problems. 

PACS numbers: 68.15.+e, 47.45.Gx, 82.39. Wj, 68.43.-h 
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Hydrophobic surfaces are at the origin of many sur- 
prising and potentially useful effects such as hydro- 
dynamic slippage at hydrophobic surfaces^, S| and the 
formation of nanobubbles at the surface [J]. A common 
feature underlying many of these phenomena is the for- 
mation of a layer of depleted water density near the sur- 
face Q - a vapor layer in the case of extremely hydropho- 
bic surfaces. Vapor-liquid interfaces have been found in 
recent spectroscopic experiments [6j and computer sim- 
ulations Q to attract large and polarizable ions such as 
bromide and iodide, but not small ions like sodium and 
chloride. This ion-specific behavior, contrary to tradi- 
tional theories of electrolyte interfaces, which only take 
into account differences in ion valency 0, is behind the 
substantial dependence on anion type of the surface ten- 
sion of aqueous solutions of halide salts |9j. Just as ion 
specificity affects equilibrium properties of vapor-liquid 
interfaces like surface tension, a similarly important role 
is expected for dynamic phenomena near the "vapor- 
liquid- like" interfaces induced by hydrophobic surfaces. 
The implications are considerable for fluid transport in 
microfluidic ("lab-on-chip") devices [10j], for which sur- 
face effects are predominant and electrokinetic techniques 
for driving flows widely used, but also for the modeling of 
biological systems for which ion-specific Hofmeister 
series are ubiquitous [11 1 . 



In this work, we investigate by computer simulations 
the anomalous electrokinetic effects that arise in electro- 
osmotic (EO) flow through hydrophobic channels due to 
interfacial ion specificity. Furthermore we develop a sim- 
ple model, comprising continuum hydrodynamic equa- 
tions and a modified Poisson-Boltzmann (PB) descrip- 
tion for the ion density distributions. Remarkably, our 
theory is able to predict accurately the effects of ion 
specificity on the simulated EO flow profiles and zeta 
potential, pointing furthermore to the crucial role of the 
hydrophobic solvation energy. Our analytic theory is a 
powerful tool for describing specific ion effects and their 



consequences on dynamics. 

The system studied comprised a solution of 
monoatomic, monovalent salt ions in water, con- 
fined between two parallel solid walls. A total of 2160 
fluid molecules were used in all cases and the SPC/E 
simple point charge model was employed for the aqueous 
solvent. Each wall was composed of 648 atoms arranged 
in three unit cell layers of an fee lattice oriented in the 
(100) direction (lateral dimensions (x,y): 48.21 A x 
32.14 A). To model charged surfaces, identical charges 
were added to each of the atoms in the top solid layer 
in contact with the fluid. The inter-wall distance was 
adjusted such that the average pressure, defined by the 
force per unit area on the solid atoms, was approximately 
10 atm in equilibrium simulations. Periodic boundary 
conditions were applied in the x and y directions, while 
empty space was added in the z direction such that the 
total system was three times as large as the primary 
simulation cell, which was centered at z = 0. 

Simulations were carried out with the LAMMPS [l2j 
molecular dynamics package. Bond length and angle 
constraints for the rigid water molecules were enforced 
with the SHAKE algorithm and a constant temperature of 
298 K was maintained with a Nose-Hoover thermostat 
(applied only to degrees of freedom in the y direction 
in the flow simulations where flow is along x). Elec- 
trostatic interactions were calculated with the particle- 
particle particle-mesh (PPPM) method, with a correc- 
tion applied to remove the dipole-dipole interactions be- 
tween periodic replicas in the z direction. Short-ranged 
van der Waals interactions between particles were mod- 
eled with the Lennard- Jones (LJ) potential, Vij(rij) — 

Asij (<Tij /r~ij) 12 — (dij /rij) 6 for an interparticle separa- 
tion of rij and particle types i and j (cr^ = (an + ajj) /2 



\f £ n £ ii)- All LJ interactions were truncated 



and e, 

and shifted to zero at 10 A. For the solid atoms, LJ pa- 
rameters were chosen to create a physically reasonable, 
albeit idealized, surface: we took a ss = 3.37 A, used a 
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FIG. 1: Simulated density profiles of negative ions (solid 
lines), positive ions (dashed lines), and water (dotted lines) 
for roughly 1-M solutions of: (a) Nal and (b) NaCl between 
neutral hydrophobic surfaces; (c) Nal at liquid-vapor inter- 
face; (d) Nal between neutral hydrophilic surfaces. 

close-packed density, p s — cr~ 3 , and chose e ss — 0.164 
and 2.08 kcal/mol to create, respectively, a hydrophobic 
and a hydrophilic surface, as characterised by the contact 
angles of a water droplet on these surfaces of roughly 140 
and 55°. 

EO flow of solutions of either Nal or NaCl were stud- 
ied, the only difference between the two cases being an- 
ion size. Except for one case, we used ion LJ parameters 
from Ref. we chose an — 6.00 A for I (instead of 
5.17 A) to reproduce approximately liquid-vapor inter- 
facial ion densities measured in simulations of Nal/ water 
solutions of similar concentration but using more com- 
plex polarizable force fields 0. Our simulated density 
profiles are shown in Fig. [1}:. Although simulations have 
shown that ion polarizability plays a role in stabilizing 
I - at the air-water interface Q, the dominant contri- 
bution to the stabilization free energy is associated with 
the solvation energy 14|, which is accounted for in our 
simulations. Thus, we regard our simple parametrization 
of the iodide LJ diameter, coupled with the use of non- 
polarizable force fields, as adequate for the purpose of 
capturing the dynamic consequences of the experimen- 
tally observed surface enhancement of I - . 

The effects of anion size and surface wettability on in- 
tcrfacial ion densities are illustrated in Fig. [TJ While 
Cl~ is not found near the hydrophobic surface in Fig.[T}D, 
Fig. UK shows a substantially enhanced interfacial T~ con- 
centration. No such enhancement is seen for I - ions near 
the hydrophilic surface (Fig. QJ1) even though the direct 
ion-solid interactions are stronger in this case, indicating 
that the ion density profiles arise largely due to the water 
structure induced by the surface. 

EO flow was induced in our simulations by applying 
an electric field E x of 0.05-0.4 V/nm in the x direction 
(linear response to the applied force was verified for all 
reported results). Starting from an initial random con- 
figuration with zero total linear momentum, simulations 
were carried out for roughly 10 ns, with statistics col- 
lected only after the steady state had been reached (typ- 
ically 1 ns). Surface charge densities £ of 0, ±0.031, and 
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FIG. 2: Top: Velocity profiles in a hydrophobic chanel for 
E = -0.062, 0, and +0.062 C/m 2 (from bottom to top) 
with (a) [Nal] « 1 M and (b) [NaCl] « 1 M. The simu- 
lation results (symbols) are compared with the resolution of 
the modified PB equation using (see text for details) the Step- 
Polarization model (dashed), and the Step-Polarization model 
with A = (dotted). Typical error bars for the theoreti- 
cal curves are shown. Error bars in the simulated velocities 
are roughly the size of the points. Bottom: ionic charge den- 
sity profile p c (z) for [Nal] « 1 M with (c) E = and (d) 
E = +0.062 C/m 2 . The symbols are from simulation. For (c) 
and (d) the lines are solutions of the modified PB equation 
with the Step-Polarization model. 

±0.062 C/m 2 and electrolyte concentrations of approxi- 
mately 0.2 and 1 M (8 and 40 ion pairs, respectively, for 
X = 0) were studied. 

The measured velocity v x (z) is shown in Fig. [5] for the 
1-M solutions in a hydrophobic chanel with E = and 
±0.062 C/m 2 ; it has been scaled by the bulk viscosity 
7?, bulk dielectric constant e, and applied electric field 
E x for ease of comparison with the zeta potential, de- 
fined in terms of the velocity in the channel center as 
£ = —rjv x (0)/(e(ieE x ), where eo is the vacuum permit- 
tivity. For e, we used the dielectric constant of pure 
SPC/E water under similar thermodynamic conditions, 
e w = 68 fl5| . The zeta potentials for all the surface 
charges are given in Fig. [3] 

Figures [2] and [3] clearly show the sensitivity of the 
EO flow to anion type, particularly for the neutral and 
positively charged surfaces. (The flow for the negatively 
charged surfaces is dominated by the excess of cations, 
Na + in all simulations) . A noteworthy point is the mea- 
surement of a non-zero £ potential (C — —40 mV) for the 
neutral hydrophobic channel with a solution of Nal, even 
though the total electrostatic force exerted on the charge- 
neutral fluid is zero. This is in strong contrast to the 
traditional theory of the electric double layer [H , though 
observed in pre vious experiments [lj| [ljj and computer 
simulations [l8j. By contrast, £ for NaCl in the same 
channel was negligible. Our measured C potentials of 
and -38 mV respectively for 1-M NaCl and Nal are con- 
sistent with experimental surface potentials of roughly 
and -20 mV respectively for vapor-liquid interfaces of the 
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FIG. 3: (a) Zeta potential of the hydrophobic surface versus 
surface charge for 1-M solutions of Nal and NaCl (simulation: 
Nal - circles; NaCl - squares). The lines are solutions of the 
Stokes equation with p c from solving PB equation using the 
Step-Polarization model (dashed) (see text) . Error bars in the 
simulated £ are roughly the size of the points, (b) Schematic 
of an ion at solid-liquid interface, illustrating the origin of 
^hyd an d its calculation. 

same solutions [19j ; £ for Nal is also of similar magnitude 
to the value of -9 mV measured by electrophoresis of neu- 
tral liposomes in 1-M KI [17] , for which ion-specific effects 
should be smaller due to the greater similarity in size of 
K + and I - compared with Na + and I - . Although not 
shown in Fig. [3j we find that the zeta potential for X = 
is not sensitive to electrolyte concentration. Also, ( was 
insignificant for Nal in the neutral hydrophilic channel. 

The anomalous result for the uncharged walls can be 
understood in terms of continuum hydrodynamics, in 
which the EO flow is described by the Stokes equation Q , 

d T z i z) = ~if PeO), where p c {z) = e [p+{z) - p_(z)] is 
the total charge density due to cations and anions of 
number density p± (z) and e is the elementary charge. 
Exploiting the symmetry of our system about z = 
and integrating the Stokes equation twice with boundary 
conditions (BCs) v x \ z=Zh = b^\ z=Zh and ^\ z=0 = 0, 
where b is the slip length applied at the hydrodynamic 
boundary Zh 



gives 



= T]V x (0) 



1 



dz'(z' - z h + b)p c (z') . (I) 



According to Eq. (jTJ) , ^ is proportional to the first mo- 
ment of the charge distribution p e relative to an origin at 
the shear plane, z s — z^ — b (the plane where the non-slip 
BC applies). Unless p e (z) = everywhere, this quan- 
tity will generally be non-zero even if the total charge, 
J z dz'p c (z'), is zero, as is the case for Nal near the un- 
charged hydrophobic wall due to the differing propensi- 
ties of Na + and f - for the surface, ft has been suggested 
that such a non-zero £ potential occurs for some non- 
charged surfaces due to ion-specific "binding" [I7| , to the 
presence of an immobile interfacial layer of charge [lij ]. 
or more generally to a reduced mobility in the interfa- 
cial layer [I^]. In contrast, our present results show that 



£ will be non-zero even if all of the charge is fully mo- 
bile. Another interesting consequence of Eq. ([T]) is that, 
as long as b is finite, surface slippage makes no contribu- 
tion to the velocity of a charge-neutral fluid containing 
only mobile charge: i.e. the system behaves as if b = 
and the flow is independent of the solid-fluid friction. As 
a matter of fact, the global fluid neutrality requires the 
wall-to-fluid force to vanish in the steady state which im- 
poses both the slip velocity and the velocity gradient to 
vanish at the wall (see fig. [2] for S = 0). In this respect, 
the case where b is infinite appears singular as the veloc- 
ity at Zh need not vanish: the velocity is determined by 
momentum conservation, as momentum cannot be trans- 
ferred to the frictionless surface, and in the end no net 
flow is achieved (not shown). 

So far, we have presented a general explanation for 
the observed ion-specific electrokinetic effects; however, 
of additional practical value would be a model capa- 
ble of quantitatively predicting the ion density and EO 
velocity profiles. With this aim, we have sought to 
construct the minimal physically accurate model for 
p e for use in the Stokes equation for v x . To obtain 
p c , we solved the one-dimensional Poisson equation 
37 [~ e ^V{z) + P(z)] = p e (z), with Neumann BCs ap- 
plied at the position z w of the surface charge and a 
mean-field approximation for the ion densities, p±{z) — 
po exp {— (3 [±eV(z) + U^ t (z)] }, where po is the bulk 
ion density and U^ t is an external potential acting 
on the ions due to interactions other than the electri- 
cal potential V. For the polarization of the medium, 
P, we assumed e(z) to display a step-function behav- 
ior at the vapor-liquid interface (from eg to e w ) so that 



P(z) = -e [e„ 



, for z > zq and P = other- 



wise, with zo the position of the first peak in the simu- 
lated water oxygen density distribution function ( "Step- 
Polarization" (SP) model). 

For the external potential, we used the sum of three 

■Ut,A- The first two 



components: Uf xt = f7^ age + U± all 



hyd 



terms, fq^ age and £/^ all , are respectively the image po- 
tential acting on the ions due to the dielectric interface 
at z Q (Eq. (3) in @) and the ion-solid LJ interaction, 
obtained by integrating the inter-particle LJ interaction 
over a uniform density p s of solid atoms occupying the 
z < z w half-plane. The final term, U^ yd , is the free en- 
ergy to create an ion-sized cavity in the fluid, i.e. to 
solvate a solute with no attraction to the solvent. This 
hydrophobic solvation energy has generally been ignored 
in calculations of interfacial ion densities, since it is negli- 
gible compared with electrostatic interactions for typical 
small ions like Na + or Cl~. We took U^ yd to be pro- 



portional to the volume V imm of the ion immersed in the 
liquid in the z > zq half-plane (see Fig.[5]3): 

U^ d (z) = Co (V± m (z) - Vt) , (2) 

with V ; ~L the total volume of the ion of solvent-excluded 
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radius r±. We took r± from bulk simulations of ions in 
water as the radius at which the ion-water radial distri- 
bution function fell to 1/e of its bulk value (2.24, 2.98, 
and 3.73 A respectively for Na + , CI", and I - ). For the 
proportionality constant, Co = 2.8 x 10 8 J/m 3 , we used 
the solvation free energy per unit volume measured un- 
der similar thermodynamic conditions in simulations of 
hard-sphere solutes of radius 0-5 A in SPC/E water [20| . 
We omitted a final plausible term in U^ t , the Born solva- 
tion energy U^ oin for charging the ion-sized cavity, as we 
found it made little difference to our results, at least us- 
ing a relatively simple expression employed by Bostrom 
et al. 0. 

The Stokes equation was solved using the calculated 
p c {z) and r\ and b measured independently in Poiseuille 
and Couette flow simulations respectively [2lj]. For E = 
0, we used b = 0, as justified above. BCs were applied in 
all cases at the position of the first peak in the sim- 
ulated water oxygen density distribution function [22TJ . 
As a test of the validity of the continuum hydrodynamic 
description, we solved the Stokes equation using the ex- 
act p e {z) from our simulations and found almost per- 
fect agreement with the simulated velocity profiles (not 
shown). Both the charge density profile p c (z) and ve- 
locity profiles calculated from the modified PB theory 
described above with the full U^ t are in good agreement 
with the simulated results, as shown in Fig. [51 The result- 
ing prediction for the £ potential reproduces very well, 
both qualitatively and quantitatively, the simulation re- 
sults, as shown in Fig. [3] It should be noted that the 
non-monotonic behavior of £ as a function of £ in Fig. [3] 
is due to the decrease in the slip length b with E. Note 
that it is possible to replace the SP model for P(z) by 
the exact value of the polarization (the gradient of which 
is equal to minus the charge density due to water in our 
simulations): doing so yields an even better agreement 
of the predicted £ potentials with simulation results (not 
shown) but at the expense of using the simulated water 
charged density profile as an input. Finally, when U^ yd is 
neglected the calculated ion density profiles and veloci- 
ties are significantly wrong for the neutral and positively 
charged surfaces (see Fig. [2]), pointing to the crucial role 
of the hydrophobic solvation energy. This is not an issue 
for the negatively charged surfaces, since the flow is dom- 
inated by Na + , for which d is negligible. Although not 
shown, we found that the conventional theory of the elec- 



tric double layer, which assumes e(z) — e w and U e 







everywhere, performed very poorly in almost all cases. 

In summary, we have shown that anomalous electroki- 
netic effects such as non-zero C potentials for uncharged 
surfaces are generic features of EO flow in hydropho- 
bic channels when the dissolved cation and anion dif- 
fer substantially in size. We have also developed a sim- 
ple model, comprising continuum hydrodynamic equa- 
tions and a modified PB description for the ion densities, 
which accurately predicts the simulated flow profiles. We 



have found that the incorporation in the model of an 
ion-size-dependent hydrophobic solvation energy, which 
favors interfacial enhancement of large ions, is crucial to 
reproducing the ion-specific effects observed in the simu- 
lations. Such an analytic theory, which is able to capture 
the subtle and complex effects of the interfacial specificity 
of ions, provides a very useful framework for the model- 
ing of biolo gica l systems, for which Hofmeister series are 
ubiquitous [ill ]. 
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